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Abstract 

A computer adapted theory for self-consistent calculations of the wavevector- and 
frequency-dependent dielectric constant for interaction site models of polar systems 
is proposed. A longitudinal component of the dielectric constant is evaluated for 
the TIP4P water model in a very wide scale of wavenumbers and frequencies using 
molecular dynamics simulations. It is shown that values for the dielectric permittiv- 
ity, calculated within the exact interaction site description, differ in a characteristic 
way from those obtained by the point dipole approximation which is usually used in 
computer experiment. It is also shown that the libration oscillations, existing in the 
shape of longitudinal time-dependent polarization fluctuations at small and inter- 
mediate wavevector values, vanish however for bigger wavenumbers. A comparison 
between the wavevector and frequency behaviour of the dielectric constant for the 
TIP4P water and the Stockmayer model is made. The static screening of external 
charges and damping of longitudinal electric excitations in water are considered as 
well. A special investigation is devoted to the time dependence of dielectric quan- 
tities in the free motion regime. 
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1 Introduction 



The study of dielectric properties of polar liquids by computer experiment is still 
a major challenge, given that the dielectric quantities are very sensitive to long-range 
intermolecular interactions and because long trajectories are required in order to obtain 
adequate statistical accuracy. For this reason, until now, the calculation of the wavevector- 
and frequency-dependent dielectric permittivity, e(k, u), in the whole region of k and uj has 
been performed for the simplest model of polar systems only, namely, for the Stockmayer 
fluid [1]. This model is related to a class of molecular models describing interactions 
by point dipoles embedded in molecules. The Stockmayer system, however, does not 
reproduce satisfactorily the dielectric behaviour of any real polar fluid. 

At the same time, more realistic interaction site (IS) models have also been consid- 
ered [2-16]. In these atomic models, intermolecular potentials are presented as a sum of 
pairwise additive site-site terms. Although IS models do not take into account internal 
degrees of freedom such as electronic polarizability and, as a rule, intramolecular vibra- 
tions, they are able to reproduce experimental results in a more satisfactory way. There 
are two approaches to investigate the dielectric properties of IS models in both computer 
experiment and theory. In the first, point dipole (PD) approach, the charged sites of 
each molecule are replaced by a point dipole, located in the molecular centre of mass, 
at constructing the reaction field and the microscopic polarization vector. The site-site 
contributions are taken into account implicitly only, namely, in the intermolecular poten- 
tials, calculating statistical averages. In the second, IS description, the reaction and the 
polarization vector are constructed in view of explicit details of charge distribution within 
the molecule [16-18]. 

Main attention in early studies was directed to evaluate the static and frequency- 
dependent dielectric constant, e(u), in the long-wavelength limit. In particular, such 
quantities have been calculated in molecular dynamics (MD) simulations for the rigid 
MCY [5] and TIP4P [6] models of water as well as for the flexible SPC potential [7] . The 
frequency dependence of the dielectric permittivity at nonzero wavevectors was investi- 
gated for the MCY [2] and TIP4P [9, 10] water models, for models of methyl cyanide 
[3, 4], methanol [11] and methanol-water mixtures [13]. But the investigations have been 
restricted to small wavevector values for which the PD approach and the IS description 
are expected to be identical. The wavevector dependence, e(k), including high wavevector 
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values, was considered in the static limit for methanol [8], methanol-water mixtures [12] 
and for the MCY water [14], using, however, the PD approximation. And only quite a 
few simulations [15, 16] have been performed to define the entire wavevector dependence 
within the exact IS description. To our knowledge, there are no computer experiment data 
concerning the entire frequency dependence of the dielectric permittivity of IS models at 
arbitrary wavenumbers. 

It is necessary to point out also the following aspect of computer simulation of IS 
models. The calculation of the dielectric constant by MD requires explicit considerations 
of a finite-size medium with periodic boundaries. Usually either the cumbersome Ewald 
summation technique [19] or the reaction field [20] is applied for treating long ranged 
interactions. As is now well established for systems of point dipoles, proper calculations 
can be made with the both methods [21]. In the case of IS models the pattern is different. 
The usual point dipole reaction field (PDRF) geometry, being exact for macroscopic 
systems of point dipoles, may not be necessarily applicably to interpret simulation results 
of IS models [22]. Recently [16], it has been shown by actual calculations that uncertainties 
of the dielectric quantities are significant if the standard PDRF geometry is used in 
computer simulations. As a result, an alternative approach has been proposed. This 
approach deals with the IS description at constructing the reaction field and leads to the 
so-called interaction site reaction field (ISRF) geometry. It has been demonstrated that 
the ISRF geometry exhibits to be much more efficient with respect to the usual PDRF 
for the investigation of dielectric properties of IS models. Within the ISRF geometry, 
the longitudinal component, £ L (k), of the wavevector-dependent dielectric constant for 
the MCY water model has been evaluated in a wide wavenumber range by Monte Carlo 
simulations [16]. Using correct microscopic variables for the polarization vector allows 
one to achieve the true infinite wavevector behaviour £ L {k — > oo) = 1. 

In the present paper we propose a self-consistent theory for the calculation of the lon- 
gitudinal wavevector- and frequency-dependent dielectric constant, £ L (k, u), for IS models 
in computer experiment. This theory is applied to the TIP4P potential using MD simula- 
tions. We show that the PD approach is valid for describing the frequency dependence of 
the dielectric constant at very small wavenumbers only. For greater wavenumber values, 
the influence of higher order multipoles becomes important. Within the IS description, 
the longitudinal permittivity is evaluated in a rather very large scale of wavenumbers and 
frequencies up to the infinite wavevector and frequency limit, where e L (k,u) — > 1. 
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2 The reaction field for IS models 



We consider a polar system with iV identical molecules composed of M interaction 
sites which are enclosed in a volume V. The microscopic field, created by the molecules at 
point r G V and time t for a macroscopic system is equal to £(r, t) = linijv- > oo Yh=i X^ii 
q a {r — rf(t)} j\r — r"(t)| 3 , where q a and rf(t) are the charge and position of site a within 
the molecule i, and the sum extends over all molecules and charged sites. This sum, 
however, can not be calculated exactly in computer experiment which deals with finite 
samples. Therefore, we must restrict ourselves to a finite set of the terms for which 
\r — rf(t)\ < R, where R is a cut-off radius. The radius R does not exceed half the cell 
length \/V if a cubic sample and the toroidal boundary conventional (tbc) are used in 
simulations. 

Now the following problem appears. How to estimate the cut-off field caused by the 
summation over the unaccessible region \r — rf(t)\ > R ? A solution of this problem can 
be found within the ISRF geometry [16]. The result for conducting boundary conditions 
is 

*(r. t) = |W - |r - r ? (t)|) ^g, - . (D 

where is the Heviside function, i.e., 0(p) = 1 if p > and 0(p) = otherwise. The 
first term in the right-hand side of (1) describes the usual Coulomb field, while the second 
contribution corresponds to the reaction field in the IS description. Performing the spatial 
Fourier transform of (1) one obtains 

£(k,t)= /dre- fc -£(r,t) = -47r(l-3^^)p L (fc,t), (2) 

V,thc 

where 



■ , N,M 



AXM) = V2 E fee-*-™ (3) 

K i,a 

is the longitudinal component for the microscopic operator of polarization density for IS 
models [18] and j^z) = —cos(z)/z + sin(,2)/,2 2 denotes the spherical Bessel function of 
first order. It is worth mentioning that as far as relativistic effects are excluded from our 
consideration, the electric field (1) does not take into account retardation corrections. For 
the same reason, we have neglected also the contribution to the field caused by dynamical 
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magnetic fields of moving charges. Because of this, the electric field has a static form and 
depends on time implicitly only. This field is pure longitudinal and can be defined uniquely 
via the longitudinal component of the polarization vector P^(k,t), that is confirmed by 
equation (2). 

Let us apply an external electric field E Q (k, t) to the system under consideration, so 
that the total field is E(k,t) = E Q (k,t) + £(k,t). The longitudinal, wavevector- and 
frequency-dependent dielectric constant is defined via the material relation 4iiP^(k, uj) = 
(e L (k,uj) - i^E L (k,u), where P L (k,u) = (P h (k,u)^ and E L (k,u) = (kk-E(k,uo)^ 
are longitudinal components of the macroscopic polarization and total field, ( ) denotes 
the statistical averaging at the presence of the external field, k = k/k, k = \k\ and 
the time Fourier transform JF(k,u) = /f^ dt e~ lujt JF(k,t) has been used for the func- 
tions Pi,(k,t) and E(k,t). The perturbation theory of first order with respect to E Q 
yields (P L (fc,w)> = -J Q dte-^-(P h (k,OyP h (-k,t)) o E (k,oj)/Vk B T, where ( ) Q 
denotes the equilibrium average at the absence of the external field, and ks and T are the 
Boltzmann's constant and the temperature of the system, respectively. Then, eliminating 
E Q (k,u) from the presented above expressions, we obtain the desired fluctuation formula 

e T (k,u)-l 9yJ?J-G L (k,t)) 

L ', > = 7 — \ v — = Qy^iJ - g T (k,t)) . 4 

e L (k,u>) l + 27y^(-G L (k,t)) Jl (kR)/(kR) V lV n 

Here 

[P L (k,oyp L (-k,t)) 

is the longitudinal component of the wavevector-dependent dynamical Kirkwood fac- 
tor for the finite sample, ii = = \J2^ Qa r i\ denotes the permanent magnitude 
of molecule's dipole moment, Jzf iw (...) = / °° ... e~ lw 'dt designates the Laplace trans- 
form, y = 4irNfi 2 jwk^T and Gi,(k,t) = dGi,(k,t)/dt. In the case R — > oo, when 
j l (kR)/(kR) — > 0, the computer adapted formula (4) reduces to the fluctuation for- 
mula for macroscopic systems in terms of the infinite-system Kirkwood factor g L (k,t) = 
Hindoo G L (k,t). 

It is essential to emphasize that the fluctuation formula (4) takes into account finite- 
ness of the system explicitly by the factor j 1 (kR)/(kR). As a consequence, for suffi- 
ciently large systems, the bulk (N, V — > oo) dielectric constant can be reproduced via 
the finite-system Kirkwood factor Gi,(k,t) calculated in simulations. However, in or- 
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5 



der to achieve this self-consistency in the evaluation of the bulk dielectric constant, the 
equilibrium averaging in Gi,(k,t) must be performed for systems with the intermolec- 
ular potential which leads exactly to the microscopic electric field £(r,t) (1) at which 
the fluctuation formula has been derived. This intermolecular potential can be deter- 
mined via the relation £(r,t) = — V$(r,t)), whence <P(r,t) = X)f^ M 0i t ( r ) t), where 
(f>t(r,t) = q a {l/p a i {t) + \p a i {tf/R z + C}, pf(t) = \r-rf(t)\ and C is, in general, an arbi- 
trary constant which for infinite systems is chosen as 0"| p «_^oo = 0. In our case, according 
to the toroidal boundary conventional, (/>f\ p a =R = and, therefore, C = —3/2R~ 1 . Then 
the intermolecular potential of interaction is q b (f>i(r b j(t)) = Z)£fb9a0j( r £(*))j i- e > 

M / i -i |„a _ r b\2 o \ 

*5 = E^e ( fl - W - (^-^ + 5 Lt_l. _ _) , (6) 

where the site-site cut-off |rf — r*| < R is performed. 

The fluctuation formula (4) has already been applied previously in the long-wavelength 
limit at the investigation of the static and frequency-dependent dielectric constant for the 
MCY [5] and TIP4P [6] models. However, acting within a semiphenomenological frame- 
work, it was not understood how to perform the truncation for intermolecular potentials. 
As a result, the molecular cut-off = |t*» — Vj\ < R, where Vi is the centre of mass of 
the ith molecule, and the usual PDRF have been suggested: 

It is easy to show that the self-consistent potential (6) can be reduced to (7) in one case 
only, namely, when d/R — > 0, where d = 2 max a \rf — T*j| is the diameter of the molecule. 
In this case, the positions for sites and the centre of mass are undistinguished within 
the same molecule. For IS finite samples, where d/R ^ 0, the PDRF potential (7) may 
affect on a true macroscopic behaviour of the system considerably. At the same time, 
the intermolecular potential (p 1 ^ corresponds completely to the conditions at which the 
fluctuation formula (4) has been obtained. Therefore, as far as this formula is applied to 
treat simulation results, the ISRF potential (6), instead of (7), must be used in computer 
experiment to reproduce a correct value for the dielectric constant. 

Nevertheless, the molecular cut-off scheme can also be acceptable, but the PDRF 
potential (7) as well as the fluctuation formula (4) need to be modified. In fact, the 
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intermolecular potential y?™ takes into account dipole contributions only into the reac- 
tion field. Additional terms can be identified within the molecular reaction field (MRF) 
approach. The result is [16]: 



0(iV j [h~^! *t ~^ ^ (<s) 

where qj = Qai^i^i ~ ^? 2 I) is the quadrupole moment of the ith. molecule with 
respect to its centre of mass, <5" = r\ — I is the unit tensor, R d = R — d/2, and 
multipoles of higher orders have been neglected. The fluctuation formula corresponding 
to the potential (8) can be derived as follows. 

The microscopic electric field created by the molecules at point r, coinciding with the 
centre of the spherical cavity of radius R, surrounded by the infinity conduction medium, 
can be presented in the MRF geometry as 

£(r, t) = £ B(R d -\r- rM) IE ga ,r~rfff,3 + ^rl ( 9 ) 



i=l \a=l 



\r-r?(t)\ 3 R 3 



(quadrupole and higher-order moments do not contribute the field in the centre of the 
cavity). Performing the spatial Fourier transform of (9) yields after some algebra 

i(k, t) = -4, { E PiH, y, t) ~ ^> ') 1 • do) 



where 

,f N M N 

k' 2 i 



P(k,k',t) = ^E^ r,(f) E? s e- ifc '- a ' w , M(k,t) =E^(*)e- ifc - ri(t) , (11) 

1 a=l i=l 



k, h' are one of the allowed wavevectors nk min of the reciprocal lattice, n designates a 
vector with integer components and k min = 2n/ \/V. As we can see P(k, k, t) = Ph(k, t), 
while M{k) is the microscopic operator of polarization density for a point dipole system. 
Its longitudinal component M L (k,t) = kk-M(k,t) coincides with P L (k,t) in the point 
dipole limit: d — > 0, q a — > oo, provided fi — > const. Using the internal field (10) and 
applying the perturbation theory with respect to the external field, similar to that as in 
the case of the ISRF geometry, we obtain the following fluctuation formula 



e L (k,u)-l 9^(-G L (M)) 



{ 6kR d R* q(q ^ 0) qR d 



, (12) 
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where the summation is extended over the infinite set nk min of nonzero wavevectors q, 
G L (k, q, t) = (P L (k, k + q, 0).P L (-fe, t))jNfi 2 , Q h (k, t) = (M L (fc, 0)-P L (-fc, t)) Q /iV / i 2 
are time correlation functions and Ph{k, k') = kk-P(k, k'). 

As was mentioned earlier, internal electric fields of classical systems are pure longitu- 
dinal. However, as it follows from the structure of equations (11), the field £(k,t) (10) 
for finite systems may contain a transverse component as well which vanishes in the limit 
R — > oo of infinite systems only. At the same time, in the ISRF geometry the field £(k, t) 
(2) remains by a longitudinal one even for finite systems. Therefore, from a physical 
point of view, the ISRF geometry is more natural than the MRF approach, because it 
does not influence on the true structure of electric fields. Moreover, within the MRF 
geometry, the fluctuation formula (12) appears to be very complicated with respect to 
the much more simple formula (4) in the ISRF approach. It requires the knowledge of 
the additional correlation functions Q L (k,t) and G L (k,q,t). While the function Q L (k,t) 
can be evaluated in the usual way, the function G^(k,q,t) at fixed t, even for spatially 
homogeneous systems, depends upon three parameters (magnitudes of k, q and the co- 
sine between them) and its calculation is a hard problem and impractical in simulations. 
Finally, the site-site cut-off truncation (6) has yet a minor advantage over the molecular 
cut-off scheme (8) because the intermolecular potential of interaction ip ls is continuous 
and continuously differentiable. This avoids the system energy drift associated with the 
passage of sites through the surface of the truncation sphere. 

3 Application to the TIP4P potential 

MD simulations have been carried out for the TIP4P potential [23] in the microcanon- 
ical ensemble at a density of 1 g/cm 3 and at a temperature of T = 293 K using the 
ISRF geometry (6). We considered N = 256 molecules in the cubic volume V = L 3 to 
which the toroidal boundary conditions have been applied. The interaction cut-off radius 
was half the cell length, R = L/2 = 9.856A. The simulations were started from a well 
equilibrated configuration for the positions of sites, obtained by Monte Carlo simulations. 
Initial velocities of the molecules were generated at random with the Maxwell distribu- 
tion. The equations of motion were integrated with a time step At = 2 fs on the basis 
of a matrix method using the Verlet algorithm in velocity form. The system was allowed 
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to achieve equilibrium for 100 000 time steps. The observation time over the equilibrium 
state was 1 000 000 At = 2 ns and each 10th time step has been considered to compute 
the equilibrium averages. In order to provide the exact conservation for the total energy 
of the system during such a rather very long observation time, the velocities of atoms were 
slightly rescaled after every 500 time steps, so that the relative total energy fluctuations 
did not exceed 0.01% over the whole run. The time correlation functions were calcu- 
lated with the time step At in the interval of lOOOAt = 2 ps. The wavevector-dependent 
quantities were investigated for k — [0, 1, ... , 300] A; min , where k mill = 2n/L = 0.319A \ 

The longitudinal component G L (k) = G L (k, 0) of the wavevector-dependent Kirkwood 
factor (5) at t — 0, obtained in the simulations within the IS description, and the corre- 
sponding function g L (k) = g L {k, 0) (4), related to the infinite system, are shown in fig. 1 
by dashed and solid curves, respectively. In the PD approximation the diameter d of the 
molecule is assumed to be sufficiently small in order to be entitled to replace the true 
microscopic operator Pi,(k, t) (3) of polarization density for an IS system by its analogue 
M L (k,t) = kY^iLik-fx i (t)e~ lk " ri( '^ for a system of point dipoles. The infinite-system 
Kirkwood factor calculated within the PD approximation, (^M L (k, 0)-M L (— k, 0)^ / iV/i 2 , 
is presented in fig. 1 by open squares. For the TIP4P molecule d = 1.837A and, there- 
fore, the PD approximation can not be used for calculating the wavevector-dependent 
dielectric quantities at k>2n/d ~ 3.4A \ Indeed, as we can see from the figure, the 
PD approach reproduces values for the Kirkwood factor satisfactorily in a small region 
of wavenumbers only, namely, at k < 2.5A \ At greater wavevectors, these values differ 
considerably from those obtained within the IS description. In the infinite wavevector 
regime {k — > oo) the Kirkwood factor tends to 1/3, when the PD approximation is used, 
whereas it approaches zero with the asymptotic behaviour g^\k) within the IS description 
(see Appendix, equation (A4)). The function g: (k) is plotted in fig. 1 by the long-short 
dashed curve. 

In order to verify explicitly that the ISRF geometry reproduces adequately the wave- 
vector dependence of dielectric quantities, we have performed also additional simulations 
for calculations of g L (k) using the generally recognized Ewald geometry and the PDRF 
potential (7). These results are included in fig. 1 as well (the open circles and dotted 
curve). The parameters i] = 5.76/ L and k max = 5 A; niin were used in the Ewald summation 
of Coulomb forces. As we can see from the figure, the ISRF method leads to results 
which are identical to those obtained within the cumbersome Ewald technique. At the 
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same time, deviations of values for the wavevector-dependent Kirkwood factor obtained 
using the PDRF potential from those evaluated within the ISRF geometry are of order 
20%. They are well exhibited at intermediate values of wavevectors. Such a situation can 
be explained by the fact that the PDRF geometry does not take into account the spatial 
distribution of charges within the molecule at constructing the reaction field and, thus, the 
precision of the calculations for wavevector-dependent dielectric quantities at k ~ 2n/d 
can not exceed d/R ~ 20%. And only for great wavevector values {k > 6 A 1 ), where 
the influence of boundary conditions is negligible, the both geometries become equivalent 
between themselves and the analytical formula (A4) can be applied here. 

To analyze the wavevector-dependent dielectric function in the low frequency limit, 
c<j — > 0, it is more convenient to rewrite the fluctuation formula (4) in terms of the infinite- 
system Kirkwood factor as follows 

(l-9yg L (k))-9yiuJ?Ug L (k,t)) 
£h[ ,U) (l-9yg L (kW + (9y)WJ??(g L (k,t)) ' 1 6) 

where the equality J£[ w (— g L (k, t)) = g L (k) —koJifiu(g L (k,t)) has been used. In particular, 
the static wavevector-dependent dielectric function e L (k) — 1/(1 — 9yg L (k)) is obtained 
from the relation (13) putting uj — 0. This function, calculated within the IS description 
and the PD approximation, is displayed in fig. 2a as circles connected by the solid curve 
and as open squares, respectively. As has been pointed out previously [6], the static 
dielectric constant of the TIP4P water at infinite wavelengths (k — > 0) appears to be 
smaller than the experimental value e Q m 80 for real water and consists about 53. This 
value has been determined on the basis of MD simulations using the PDRF geometry. In 
our calculations within the exact ISRF geometry we have found that e Q = e h (k — > 0) ss 50. 
Therefore, the TIP4P water reproduces the static dielectric constant of real water even 
somewhat worse than this has been established earlier. The function e L (k) in the PD 
approximation behaves like that for the Stockmayer fluid [21, 24]. For example, in the 
infinite wavevector limit, the dielectric constant tends to the wrong Onsager limit value 
1/(1 - S y ) = -0.0649. 

The true wavevector behaviour for the dielectric constant can be reproduced within 
the IS description only. As we can see, the static dielectric constant of the TIP4P water 
is negative over a wavevector range bounded by two singularities, where 9yg L (k) — > 1 and 
e L (k) — > ±oo. The first singularity is achieved at k — k T ~ 0.297A \ whereas the second 
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singularity may be estimated [14] via the relation 



k 



AnN m 

n~^^-\Vk^h q «- (14) 



This relation is derived solving the equation 9ygf\k) = 1 and neglecting intersite contri- 
butions (the second term in the right-hand side of (A4)) in the intramolecular part gf\k) 
of the Kirkwood factor. Applying the relation (14) yields for the TIP4P molecule k n ~ 
19.71A" , that is in an excellent accord with the simulation result 20.28A . At k > k T , 
the function £ L (k), remaining positive, slightly approaches unity at k <: 100 A 1 . In the low 
frequency limit at k* = k l , k n we put 1— 9yg L (k*) = and obtain from (13) that the dielec- 
tric constant behaves as e L (k*,oj) = —i/ujTl OI (k*), where rl° r (k) = / °° g h (k,t)/g L (k)dt is 
the correlation time. Therefore, the real part of the dielectric permittivity in this regime 
is equal to zero, whereas the imaginary part tends to infinity as 1/uj. This result can be 
better understood introducing the generalized coefficient a L (k,cu) of polarization conduc- 
tivity, which defines the macroscopic current Ii,(k,t) = dPi,(k,t)/dt in the frequency 
representation as I^(k,u) = a L (k,ou)E^(k,u}). Taking into account that in this represen- 
tation Ih(k, uj) = iuP^(k, uj), we can express the polarization conductivity in terms of the 
dielectric constant as a L (k,u) = f^(e L (k,uj) — 1). Thus in the limit of small frequencies 
we have that cr h (k) = lim a L (k, uj) = lj AixT™ T (k) ^ if k G k* and a h (k) = otherwise. 

In view of the existence of a nonvanishing coefficient for the generalized conductivity 
in the static limit, the following question arises. Does it violate the well-known law that 
static macroscopic currents are absent in the dielectrics? We can be sure that this law re- 
mains in force, because & L {k) ^ when |e L (A;)| — > oo and, therefore, the total electric field 
in the system vanishes because of E L (k) = kk-E Q (k) /s L (k) — > 0, so that the polariza- 
tion current does not appear. In this case, the longitudinal external field is compensated 
completely by the internal field of polar molecules. Moreover, the singularities in the 
dielectric permittivity do not lead to singularities in physically observing quantities. This 
is so because these quantities are expressed through the external electric field using multi- 
pliers of type l/e L (k) or (£ L (k) — l)/e L (k) which are free of singularities. For instance the 
ratio S(k) j (kk-E Q (k)) of the macroscopic static field of polar molecules and the external 
field is determined by the factor —(e L (k) — l)/e L (k) = —9yg L (k). The sign minus shows 
that the field of molecules is opposite to the external field, because molecular dipole mo- 
ments align always along external fields. The maximum magnitude of this factor ranges 
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up to about 44.9 at k m 3.03A (the first maximum of the Kirkwood factor), where the 
dielectric constant takes its almost maximum value —0.0228 within the negative region. 
The value 3.03A 1 corresponds to the wavelength A = 2ir/k 2.07A ~ d. Therefore, 
optimal conditions (maximal average torques acting on molecules due to electric fields) 
for alignment of dipoles along external fields are observed when spatial inhomogeneity of 
these fields has a characteristic length of its varying in coordinate space of order of the 
molecular diameter (the so called spatial resonance). A similar wavevector behaviour of 
e L is exhibited for a central force model of water [15]. 

The knowledge of the dielectric constant in the whole region of wavenumbers allows 
one to solve the problem of static screening of external charges in water. An external 
charge, enclosed in the dielectric, causes a polarization of the system. The potential of 
the electric field, created in such a way, can be written in the ^-representation as (p(k) = 
ip( \k)/e L (k), where ip(°\k) = An/k 2 is the Fourier transform for the potential (p(°\r) = 
1 jr of an unit charge in the vacuum. Then applying the inverse Fourier transform to the 
function <p(k) one obtains <p(r) = ^ j V (k)e ik ' r dk = f / ^y^dA; = ip(°\r) /e L (r). 
The function l/e L (r), which is plotted in fig. 2b, describes the static screening. It can 
be shown easily, using the asymptotic behaviour (A4), that the dielectric constant in the 
infinite wavenumber regime is expanded over inverse wavevectors as e L (k) = 1 + 1/ {r D k) 2 + 
0(k~ 4 ), where r D = l/k n . The expression e L (k) = 1 + l/(r D /c) 2 is often used in the whole 
range of wavevectors, considering processes of static screening in plasma. This leads to the 
well-known Debye exponential screening l/e L (r) = exp(— r/r D ) with the Debye's radius 
r D . In the case of dielectrics the pattern is different. The function l/e L (r), starting 
from 1 at r = 0, exhibits a pronounced oscillatory feature, reflecting the influence of the 
microscopic structure of the system. And only beginning from distances to the external 
charge of order r > 15A this function tends to its value l/e in the macroscopic limit. 

Examples for the normalized autocorrelation functions $ L (/c,t) = G\,{k,t)/G]_,{k) (5) 
and (f> L (k,t) = g L (k,t)/g L (k) (4), describing the time decay of longitudinal polarization 
fluctuations in the finite (N = 256) and infinite systems are plotted in fig. 3 by the cir- 
cles and solid curves, respectively (the infinite-system functions are obtained applying 
the inverse Laplace transform to the relation (4)). The curves are for a set of wavevec- 
tors accessible in the simulations. For the purpose of comparison, analogous functions, 
(M L (k,0)-M L (-k,t)) /(M L (k,0)-M L (-k,0)) , obtained within the PD approxima- 
tion for the finite system, are also presented in this figure by dashed curves. The self 
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(i = j in equation (Al)) part 4> (s '(k,t) = g^ s ' (k,t)/ g (k) of T (k,t) and the normalized 
function <f)f\k,t) = gf\k , t) / gf\k) (note that gl s \k) = gf\k) in the static limit) cor- 
responding to a noninteracting system, calculated with the help of stochastic simulations 
according to equation (A2), are included in fig. 3 as well and shown by the dotted and 
long-short dashed curves, respectively. As we can see, the self portion of g L (k,t) begins 
to dominate already at k > 10A; niin and the free motion regime starts at k > 40k m m- 
The normalized function g^\k,t), evaluated via the analytical formula (A3), is plotted 
in insets (e)-(h) of fig. 3 by open circles. 

The time correlation functions obtained in the IS description are identical to those 
within the PD approximation in the zero wavevector limit only but they differ from one 
another significantly at greater (k > fc m i n ) wavevector values. For example, in the infinite 
wavevector limit (see fig. 3g,h), the PD functions exhibit the pure Gaussian behaviour 
f G (k,t), whereas the IS functions have a more complicated structure of the time depen- 
dence (A3). Taking into account the remarks concerning the behaviour of these functions 
in the static limit (see fig. 1), we conclude that the PD approximation, being exact for 
systems of point dipoles, is unsuitable, in general, to investigate the wavevector- and 
frequency-dependent dielectric constant of IS models of polar fluids. 

The most striking feature in the time behaviour for the longitudinal component of the 
infinite-system function <j) L (k,t) is damped oscillations superimposed on the exponential. 
They are called librational oscillations and can be observed for the self part 4>^ s \k,t) as 
well. To show these oscillations in more detail, we have also presented 3D plot of the 
function <f> L (k,t) in fig. 4. The librational oscillations appear for systems with a slowly 
relaxing character of intermolecular torques and describe the rapid motion of a molecule 

in the average electric field ~ r z e {r) °^ ^ s neighbours (see fig. 2b). It is not a surprise, 

L ff) 
the librations are absent in the shape of the free- motion functions (j) (k,t) which can be 

related to polarization fluctuations of a dilute polar fluid. The librational oscillations have 

been found previously for the MCY [2, 5] and TIP4P potentials [6] and for a model of 

methyl cyanide [3] . However, the investigation of the librational oscillations was restricted 

to zero and very small values (~ k min ) of wavevector. In our study, performed in a rather 

very large scale of wavenumbers, we have identified the librational oscillations at small 

and intermediate wavenumber values, namely, at k < 10/c m i n ~ 3A 1 (see figs. 3, 4). They 

vanish at bigger wavevectors when the time correlation functions behave like those for a 

noninteracting system. 
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It is necessary to underline that all the functions after a sufficiently long period begin 
to decay purely exponentially in time as ~ exp(— t/r£ el (/c)), where r£ el (/c) is the relaxation 
time. The wavevector-dependent relaxation r£ el (/c) as well as the correlation r£ or (/c) times 
for the finite and infinite systems are shown in fig. 5. At k = we have found for the 
finite-system function r rel (0) = 6.7ps and r cor (0) = 6.4ps. As was pointed out above, 
the correlation functions have been evaluated in the finite time interval t e [0, 2ps]. 
The relaxation times will be used by us to calculate contributions at time integration for 
t > 2ps, where all the correlation functions exhibit almost purely relaxation behaviour. 
We indicate the longest tails in the time correlation functions for the infinite system at 
k ~ 1.7 and 2.9A \ These functions remain positive anywhere in time space, contrary to 
the Stockmayer fluid with the dipolaron behaviour [1, 25] of dipole moment fluctuations. 

Now, we are in a position to discuss the result of calculations (4) for the longitudinal 
wavevector- and frequency-dependent dielectric constant e L (k,u>) of the TIP4P water. 
Real z' L (k, u>) and imaginary s'^(k, cu) parts of £ L (k, u) = z' L (k, oj) — is'^(k, u) are shown in 
figs. 6 and 7 as solid and dashed curves, respectively. The calculation of the frequency- 
dependent dielectric constant for the TIP4P water has already been dealt with in the paper 
[6]. In this paper, however, the dielectric relaxation was investigated at zero wavevector 
value only. Moreover, the simulation results here were actually deduced within the PDRF 
potential rather than the exact ISRF geometry. As has been concluded previously [16], 
investigating the wavevector dependence of the dielectric constant of the MCY model, 
and shown above for the wavevector-dependent Kirkwood factor g L (k) of the TIP4P 
water (see fig.l), it is risky to use the PDRF geometry for systems of hundreds molecules. 
Although it is hard to tell except by numerical computations whether the results [6] for the 
frequency dependence of the dielectric constant differ significantly from the exact results 
corresponding to the macroscopic system. For this reason we have repeated the calculation 
of the frequency- dependent dielectric constant e{uS) = lim^ e L (k, uj). The result of these 
calculations, performed in the ISRF geometry, is shown in fig. 6a. Comparing the previous 
result (cf. fig. 6 of ref. 6) with our one, we can observe that apart from the slight differences 
in the low frequency regime the agreement is quite good for higher frequencies. Therefore, 
the dielectric constant is less sensitive to boundary effects at intermediate and high values 
of frequency. 

Since all the correlation functions decay exponentially at great times, in the low fre- 
quency regime the dielectric permittivity behaves like the Debye dielectric for arbitrary 
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wavevectors. Deviations can be visible at higher frequencies and at small and interme- 
diate wavevector values, where the collective molecular librations take an important role 
in forming the polarization fluctuations. Above some ou—1-10 THz, where THz=ps~ 1 , 
the relaxation process gradually changes into a resonance process characterized by an fre- 
quency of 60-200 THz (which depends on wavevector) and reflecting the rapid librational 
motion of the molecules. The main differences at zero wavevector between the frequency- 
dependent dielectric constant for the TIP4P water, Debye and Stockmayer models as 
well as real water have been done already [6, 26]. Now we consider differences in the 
behaviour on frequency between the longitudinal dielectric constant of the TIP4P water 
and the Stockmayer model [1] at nonzero wavevector values. 

Despite a similar overall shape of e L (k,cu) to that of the Stockmayer fluid, additional 
features emerge for the TIP4P water, namely, a second maximum in the imaginary part 
e''{k,uj) at k < 7k niin and a more complicated structure for the real part e' (k,u) at 
k < 20A; min . For example, the real part £' L (k, oj) for the Stockmayer system always increases 
at increasing frequency, while we can observe clearly additional maximum and minimum 
for the TIP4P water in the librational resonance range. Although the TIP4P model does 
not reproduce the dielectric properties of water at all well, since it does not account for 
polarizability and the coupling of inter- and intramolecular motions, a similar structure 
in the dielectric constant of real water for nonzero wavevector values may be present 
as well, because in real water the librational frequency is well separated from the lowest 
intramolecular vibration [27]. It is interesting also to point out that the dielectric constant 
of the TIP4P water at frequencies oj > 10 THz practically does not depend on wavevector 
in the range < k < 2k niin and even up to k = 10/c min at oj > 200 THz and remains 
the same as in the low wavevector limit (cf. fig. 6b-e). The wavevector- and frequency- 
dependent dielectric constant at great wavenumber values is shown in fig. 7. As we can 
see, the correct value of nonpolarizable systems in the infinite wavevector and frequency 
limit, = 1, can only be obtained if the dielectric constant is known up to about 
k ~ 150A; min « 50A 1 and oj ~ 1000 THz. It is worth to remark also that at k > 20/c min 
the frequency dependence of the dielectric constant can be reproduced quantitatively 
using the analytical formula (A3) for the time correlation functions of a noninteracting 
system. The dielectric constant, calculated in such a way, is shown in fig. 6f-h and fig. 7 
by circles. 

The dielectric constant e L (k,oj) at fixed real values of wavevector and frequency de- 
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scribes the response of the system on longitudinal electric fields in the form of monochro- 
matic plane waves, i.e., when E(r,t) ~ e i(wt-fe-r)_ Arbitrary inhomogeneous fields can 
be presented by a set of the plane waves via the time and spatial Fourier transform. 
Therefore, the response on such fields can be determined integrating the correspond- 
ing monochromatic contributions over the frequency and wavenumber spectrum using 
e L (k,u). The necessity of introducing pure imaginary frequencies uj = iu* with uj* > 
arises in the problem of finding electric fields E(r,t) ~ e (-w*t-ifc-r) - m ^ e system when 
the external field is turned off. Such damping in time electric fields can be obtained 
using the dispersion relation e L (k,uj) = [28]. It can be seen from figs. 6, 7 that 
this relation is not satisfied at real frequencies. So that, contrary to transverse electro- 
magnetic waves, longitudinal monochromatic waves can not propagate without external 
sources of charges. From the fluctuation formula (4) it follows that £ L (k, uj) — > when 
=^L( — 9 ' L (k,t)) — > oo. The last limit can be achieved at imaginary frequencies uj = iu* 
only, where uj*(k) — > l/rl cl (k). Thus, the electric fields are damped exponentially with 
the characteristic interval of order of the relaxation time. Maximal life times of the lon- 
gitudinal electric excitations are observed at spatial inhomogeneity of k ~ 1.7A 1 (see 
fig. 5) and consists about 1.6ps. Strong spatially inhomogeneous fields (k — > oo) disappear 
immediately after disappearing external supporting fields. 

Finally, we consider some aspects of applying the well-known Kramers-Kronig relations 
to the longitudinal wavevector- and frequency-dependent dielectric constant. In general, 
these relations connect real and imaginary parts of the function £(a>) = — 
i£"(u;) of frequency as 

7T J X 2 — UJ 1 TV J X — UJ 



From the mathematical point of view the relations (15) are fulfilled if the function 
can be presented in the form Jz? ia; (£(£)), where £(t) is a real analytical function of time. 
Then, taking into account the fluctuation formulas (4), it can be shown that in two cases 
at least, namely, when = £ L (0, to) — 1 = e(lu) — 1 and = (s L (k, uj) — j £ L (k, uj) 
for arbitrary k, the Kramers-Kronig relations will be satisfied. At the same time it is not 
obvious that these relations must be valid for the function £(a>) = e L (k,oj) — 1 at k ^ 0. 
We have established on the basis of our numerical analysis that the Kramers-Kronig 
relations can apply to the longitudinal dielectric constant in the form e L (k,u>) — 1 at 
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such wavenumbers only, where the static dielectric constant e L (k) takes positive values, 
i.e., when < k < k l and k > k lY Physically, the impossibility to use the function 
e L (k,uj) — 1 for the transformations (15) at arbitrary wavenumbers lies in the fact that, in 
fact, the external field E Q , but not the total field E, can be considered as an independent 
parameter which causes time varying all observing quantities in the system. 

4 Conclusion 

A fluctuation formula suitable for the self-consistent calculations of the wavevector- 
and frequency-dependent dielectric constant for interaction site models of polar systems 
has been rigorously derived in the situation that is typical for computer simulations. Using 
this formula, the longitudinal component of the dielectric constant has been evaluated by 
the MD method for the TIP4P model of water in a very wide scale of wavenumbers and 
frequencies. The most striking feature of interaction site models consists in existing of the 
libration oscillations in the shape of longitudinal time-dependent polarization fluctuations 
and this feature is reproduced by the TIP4P potential as well. We have showed, however, 
that the libration oscillations vanish at increasing wavenumber. Choosing the correct 
microscopic variable for polarization density, which corresponds completely to interaction 
site models, allows us to investigate the frequency dependence of the dielectric constant 
for arbitrary wavevector values. At the same time, it has been corroborated by the 
explicit calculations that the PD approximation is unsuitable, in general, for evaluating 
the wavevector- and frequency-dependent dielectric constant of interaction site models of 
polar fluids. 

Since it has now been shown that the calculation of the wavevector- and frequency- 
dependent dielectric constant in computer simulations for a given interaction site model 
is practical in principle, we believe that this fact will stimulate further research in both 
theory and pure experiment. A next paper of this series will be devoted to a more 
complicated case, namely, to the calculation of the transverse wavevector- and frequency- 
dependent dielectric constant for the TIP4P water. 

The author would like to acknowledge financial support of the President of Ukraine. 
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Appendix 



We now consider the time dependence of the longitudinal wavevector-dependent Kirkwood 
factor 

-, / N;M v 

<t<M)-^(E<^ef*™^<»>) (Al) 

in the free motion regime. In this case, molecules are statistically independent of one an- 
other and, therefore, nonzero contributions to g L (k,t) give only the self (intramolecular) part 
of terms with coincident molecular indexes (i = j) of summation (Al). The site velocities of 
molecules can be presented as vf(t) = Vj + f2i(t)x5f(t), where Vj and i?« are the transla- 
tional and angular velocities of the ith molecule, respectively. We note that according to the 
Euler equations, the angular velocities depend on time even for free rotational motion, so that 
rf(t) = n (0) + Vi t+Sf(t), where Sf(t) = 6? (0)+/ *[i7 i (t')x5f(t')]dt'. Then, taking into account 
that the translational and angular velocities of each molecule are distributed independently at 
equilibrium, the expression (Al) transforms into 

, N 1 M \ 

s 'l <*•«) = £ (e"" v ''}„ E «^ W, -' !M ■ (A2) 

" i=l \ a,b I q 

The first averaging in (A2) can easily be evaluated applying the Maxwell distribution 
fo(Va) = ( k 2tv^t ) ^ ex P (" ^t ) w ^* n res P ec t to the cartesian components (a = x,y,z) of the 
translational velocity V, where m is the mass of the molecule. As a result, we obtain ( elfc V *) = 

/ / I-oo dy - d W4,(^)/ (W (^)^ = / G (M)- 

The second averaging in (A2) is not reduced, in general, to an analytic form. However, due to the 

fact that the Gaussian multiplier f Q (k,t) quickly decays to zero with increasing t at k / 0, we 
can cast the solution for 6"(i) into the Taylor series by writing 5"(i) = 5" + [12jX($f]i + 0(t 2 ), 
where i?j and <5" are taken at time t = and terms of order t 2 and higher have been ne- 
glected. The integration over angular velocities it is convenient to perform in the moving coor- 
dinate system XYZ, attached to the molecule, in which the Maxwell distribution has the form 
/ (^o) = ( 2wk B T ) ex P ( ~ 2k B T ) ' wriere J a denote the moments of inertia along principal axes 
a = X, Y, Z. Finally, performing averaging over orientations of the molecule with respect to the 
laboratory frame xyz, we obtain the desired result 

^ (M) =^^E^ 6 / 7r sinM0 l^cos^cUcosSOexp (-^k 2 t 2 

^ a,b ® " 

^ f (A y cos - A| sin 9 sin (f>) 2 + (A^ cos 6 - A| sin 9 cos <f>) 2 
+ ( A^ sin 6» sin ^ - A Y sin (9 cos <^>) 2 1\ + ^ / t 2 \ 
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where (A%-,Ay, A^) designate the principal components of Sf and d a b = \Sf — are the 
distances between charges a and b within the same molecule. 

At a given molecular geometry the two-dimensional integral (A3) can easily be calculated 
numerically. In the static limit (t = 0) this integral is taken analytically, 

(f), n 1 1^ 2 , r sm(kd ab )\ 

It is worth to remark that the formulas (A3), (A4) at sufficiently large values of wavevector k 
can be applied to interacting systems as well, because then both intermolecular terms in (Al) 
and terms in (A3) caused by interactions (nonlinear in t) are small. 
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Figure captions 



Fig. 1. The longitudinal wavevector-dependent Kirkwood factor of the TIP4P water. The 
result obtained for the finite and infinite systems within the IS description using the ISRF 
geometry is plotted by dashed and solid curves, respectively. The corresponding result of the 
point dipole approximation for the infinite system is presented by open squares. The infinite- 
system Kirkwood factor, evaluated in the Ewald and PDRF geometries, is shown as open circles 
and dotted curve, respectively. The long-short dashed curve corresponds to the Kirkwood factor 
in the free motion regime (equation (A4)). 

Fig. 2. The longitudinal wavevector-dependent dielectric constant (a) of the TIP4P water 
within the IS description (circles connected by the solid curve) and the PD approximation (open 
squares). Two vertical lines indicate the positions of singularities. The function, describing the 
static screening of external charges in water, is shown in (b). 

Fig. 3. The normalized, time autocorrelation functions of the longitudinal polarization 
fluctuations in the TIP4P water for the finite (circles) and infinite (solid curves) systems in the IS 
description as well as for the finite system in the PD approximation (dashed curves). The dotted 
curves show the self parts of the IS functions. The correlation functions of a noninteracting 
system, calculated using the exact relation (A2) and the short time approximation (A3), are 
shown by the long-short dashed curves and open circles, respectively. We note that five of the 
six dependencies are undistinguished in (h). 

Fig. 4. The 3D plot of the wavevector- and time-dependent longitudinal polarization fluc- 
tuations in the IS description for the infinite system. 

Fig. 5. The wavevector-dependent relaxation (circles) and correlation (squares) times for 
the longitudinal polarization fluctuations of the finite system. The results for the infinite systems 
are shown by the corresponding solid curves. 

Fig. 6. The longitudinal wavevector- and frequency-dependent dielectric constant of the 
TIP4P water. The real and imaginary parts are shown as solid and dashed curves, respectively. 
Note the logarithmic scale in (a). Open squares connected by the long dashed curve in (b)-(e) 
reproduce the real part of the frequency-dependent dielectric constant at zero wavevector. The 
results obtained in the free motion regime (equation (A3)) are plotted in (f)-(h) as full (real 
part) and open (imaginary part) circles, respectively. 

Fig. 7. The longitudinal wavevector- and frequency-dependent dielectric constant of the 
TIP4P water at great wavenumber values. Other features as for fig. 6. 
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